Table of Contents
You have already seen the basics of the general linear model (linear regression, t-tests, and ANOVAs). Now, we are going to get into some expanded topics related to the general linear model: centering, interactions, and mediation.
Main Ideas
Centering – Makes interpretable intercepts
Interactions – Adds context to models
Mediation – Adds paths to models
Centering
Think back to what the intercept means in the context of a linear regression model: it is the value of the DV/outcome when everything else is at 0. There are times where that makes sense (you can have 0 dollars, but I don’t recommend it).
What if we want to have a meaningful intercept?
library(data.table)
crimeScore <- fread("http://nd.edu/~sberry5/data/crimeScore.csv")
uncentered <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT, data = crimeScore)
summary(uncentered)
Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT, data = crimeScore)
Residuals:
Min 1Q Median 3Q Max
-255.974 -44.233 8.767 42.285 223.767
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 275.4925 0.3175 867.70 < 2e-16 ***
NARCOTICS_ARR_CNT 0.7408 0.1204 6.15 7.76e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 60.39 on 92805 degrees of freedom
(305877 observations deleted due to missingness)
Multiple R-squared: 0.0004074, Adjusted R-squared: 0.0003967
F-statistic: 37.83 on 1 and 92805 DF, p-value: 7.758e-10
plot(uncentered$fitted.values, uncentered$model$SSL_SCORE)
plot(uncentered$model$SSL_SCORE, uncentered$model$NARCOTICS_ARR_CNT)
Now we can do some very simple centering: we can just subtract the mean of the item from every observation of that item.
crimeScore[, narcCenter := NARCOTICS_ARR_CNT - mean(NARCOTICS_ARR_CNT,
na.rm = TRUE)]
# Notice how I didn't assign crimeScore back to crimeScore?
# The skull operator handled it for us!
centeredMod <- lm(SSL_SCORE ~ narcCenter,
data = crimeScore)
summary(centeredMod)
Call:
lm(formula = SSL_SCORE ~ narcCenter, data = crimeScore)
Residuals:
Min 1Q Median 3Q Max
-255.974 -44.233 8.767 42.285 223.767
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 277.0178 0.1982 1397.34 < 2e-16 ***
narcCenter 0.7408 0.1204 6.15 7.76e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 60.39 on 92805 degrees of freedom
(305877 observations deleted due to missingness)
Multiple R-squared: 0.0004074, Adjusted R-squared: 0.0003967
F-statistic: 37.83 on 1 and 92805 DF, p-value: 7.758e-10
plot(centeredMod$fitted.values, centeredMod$model$SSL_SCORE)
plot(centeredMod$model$SSL_SCORE, centeredMod$model$narcCenter)
We can see that nothing really changes except for the intercept’s coefficient. Let’s think about what is being said here: with narcotics arrest count being a 0, we would expect the SSL score to be 277 on average. Since we centered that narcotics variable, though, the 0 is actually the mean of that variable (2.06).
Standardizing
We will get back to more interesting data in a second, but let’s look at some mtcars data first (we can’t possibly learn anything new about it):
cor(mtcars[, c("hp", "disp", "wt")])
hp disp wt
hp 1.0000000 0.7909486 0.6587479
disp 0.7909486 1.0000000 0.8879799
wt 0.6587479 0.8879799 1.0000000That is a pretty strong correlation. Let’s see how they behave in a model:
testMod <- lm(mpg ~ hp + disp + wt, data = mtcars)Since those variable had pretty strong correlations, let’s look at the variance inflation factor for each:
\[\frac{1}{1-R^2_i}\]
summary(lm(hp ~ disp + wt, data = mtcars))
Call:
lm(formula = hp ~ disp + wt, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-53.830 -33.345 -3.866 15.445 155.541
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.8443 31.8023 2.165 0.038776 *
disp 0.5388 0.1350 3.990 0.000411 ***
wt -14.4453 17.1038 -0.845 0.405268
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 42.85 on 29 degrees of freedom
Multiple R-squared: 0.6346, Adjusted R-squared: 0.6094
F-statistic: 25.18 on 2 and 29 DF, p-value: 4.575e-07We can plug our \(R^2\) value in:
1 / (1 - .6346)
[1] 2.736727VIF starts at 1 and goes up from there. Anything over 5 would indicate worrisome multicolinearity.
car::vif(testMod)
hp disp wt
2.736633 7.324517 4.844618 When modeling these main effects, nothing is going to change this problem. If we are going to be using our variable as interactions, we might want to standardize these variables – we would center our items and then divide by the standard deviation of the item:
mtcarsDT <- as.data.table(mtcars)
scaleVars <- c("hp", "disp", "wt")
mtcarsDT[,
(scaleVars) := lapply(.SD, function(x) scale(x)),
.SDcols = scaleVars]
summary(mtcarsDT)
mpg cyl disp hp
Min. :10.40 Min. :4.000 Min. :-1.2879 Min. :-1.3810
1st Qu.:15.43 1st Qu.:4.000 1st Qu.:-0.8867 1st Qu.:-0.7320
Median :19.20 Median :6.000 Median :-0.2777 Median :-0.3455
Mean :20.09 Mean :6.188 Mean : 0.0000 Mean : 0.0000
3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.: 0.7688 3rd Qu.: 0.4859
Max. :33.90 Max. :8.000 Max. : 1.9468 Max. : 2.7466
drat wt qsec vs
Min. :2.760 Min. :-1.7418 Min. :14.50 Min. :0.0000
1st Qu.:3.080 1st Qu.:-0.6500 1st Qu.:16.89 1st Qu.:0.0000
Median :3.695 Median : 0.1101 Median :17.71 Median :0.0000
Mean :3.597 Mean : 0.0000 Mean :17.85 Mean :0.4375
3rd Qu.:3.920 3rd Qu.: 0.4014 3rd Qu.:18.90 3rd Qu.:1.0000
Max. :4.930 Max. : 2.2553 Max. :22.90 Max. :1.0000
am gear carb
Min. :0.0000 Min. :3.000 Min. :1.000
1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
Median :0.0000 Median :4.000 Median :2.000
Mean :0.4062 Mean :3.688 Mean :2.812
3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
Max. :1.0000 Max. :5.000 Max. :8.000
scaledMod <- lm(mpg ~ hp + disp + wt, data = mtcarsDT)
summary(scaledMod)
Call:
lm(formula = mpg ~ hp + disp + wt, data = mtcarsDT)
Residuals:
Min 1Q Median 3Q Max
-3.891 -1.640 -0.172 1.061 5.861
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.0906 0.4665 43.067 < 2e-16 ***
hp -2.1362 0.7841 -2.724 0.01097 *
disp -0.1161 1.2827 -0.091 0.92851
wt -3.7190 1.0432 -3.565 0.00133 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.639 on 28 degrees of freedom
Multiple R-squared: 0.8268, Adjusted R-squared: 0.8083
F-statistic: 44.57 on 3 and 28 DF, p-value: 8.65e-11Interactions
Sticking with some boring data for just a little while longer:
library(ggplot2)
ggplot(mtcars, aes(mpg, hp)) +
geom_point() +
geom_smooth(method = "lm") +
theme_minimal()
We can clearly see the effect of horsepower on mpg.
Let’s add some additional context to this visualization:
ggplot(mtcars, aes(mpg, hp, color = as.factor(cyl))) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
theme_minimal()
We have already been doing stuff like this, but not really thinking about the model that could test it.
Let’s start by looking at a multiple regression with 2 predictor variables:
twoVars <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT + NARCOTICS_ARR_CNT,
data = crimeScore)
summary(twoVars)
Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT + NARCOTICS_ARR_CNT,
data = crimeScore)
Residuals:
Min 1Q Median 3Q Max
-194.221 -28.029 -2.221 26.550 187.586
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 291.5470 1.5197 191.850 <2e-16 ***
WEAPONS_ARR_CNT 18.0598 1.0118 17.849 <2e-16 ***
NARCOTICS_ARR_CNT 2.8072 0.2933 9.573 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 49.69 on 9037 degrees of freedom
(389644 observations deleted due to missingness)
Multiple R-squared: 0.04265, Adjusted R-squared: 0.04243
F-statistic: 201.3 on 2 and 9037 DF, p-value: < 2.2e-16Let’s explore interactions (moderation to some). A key idea here is that interactions change the nature of our model. Instead of supposing that the predictor variables act in isolation on the DV, the interaction is essentially providing expanded context for our model. We are essentially saying that the values of one variable will have an effect with values of another variable on the DV. Here is what this looks like in words:
“What is the effect of X on Y?”
“It depends on Z’s value.”
intMod <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT * NARCOTICS_ARR_CNT, data = crimeScore)
summary(intMod)
Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT * NARCOTICS_ARR_CNT,
data = crimeScore)
Residuals:
Min 1Q Median 3Q Max
-194.231 -27.995 -2.231 26.532 187.532
Coefficients:
Estimate Std. Error t value
(Intercept) 292.1106 2.2340 130.754
WEAPONS_ARR_CNT 17.5941 1.6896 10.413
NARCOTICS_ARR_CNT 2.5461 0.8134 3.130
WEAPONS_ARR_CNT:NARCOTICS_ARR_CNT 0.2173 0.6312 0.344
Pr(>|t|)
(Intercept) < 2e-16 ***
WEAPONS_ARR_CNT < 2e-16 ***
NARCOTICS_ARR_CNT 0.00175 **
WEAPONS_ARR_CNT:NARCOTICS_ARR_CNT 0.73069
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 49.69 on 9036 degrees of freedom
(389644 observations deleted due to missingness)
Multiple R-squared: 0.04266, Adjusted R-squared: 0.04234
F-statistic: 134.2 on 3 and 9036 DF, p-value: < 2.2e-16This is the model that we have estimated:
\[score = b_0 + b_1(weapons) + b_2(narcotics) + b_3(weapons*narcotics)\]
The interpretation of our main effects don’t really change.
Our interaction terms (\(b_3\)) is providing the amount of change in the slope of the regression of score on weapons when narcotics changes by one unit. So as narcotics arrests increase, we see an increase in the effect of weapons arrests on score (but only a tiny one).
To predict what the score value would be for certain values of weapons arrests, we could reformulate our model as:
\[score = b_0 + (b_1 + b3*narcotics)weapons + b_2(narcotics)\]
Before we look at those interactions, let’s compare how our models fit:
anova(twoVars, intMod)
Analysis of Variance Table
Model 1: SSL_SCORE ~ WEAPONS_ARR_CNT + NARCOTICS_ARR_CNT
Model 2: SSL_SCORE ~ WEAPONS_ARR_CNT * NARCOTICS_ARR_CNT
Res.Df RSS Df Sum of Sq F Pr(>F)
1 9037 22311153
2 9036 22310861 1 292.56 0.1185 0.7307We can compare some information about the residuals to compute an F statistic and the accompanying p-value – these models are no different from each other.
stargazer::stargazer(twoVars, intMod, type = "html", header = FALSE)| Dependent variable: | ||
| SSL_SCORE | ||
| (1) | (2) | |
| WEAPONS_ARR_CNT | 18.060*** | 17.594*** |
| (1.012) | (1.690) | |
| NARCOTICS_ARR_CNT | 2.807*** | 2.546*** |
| (0.293) | (0.813) | |
| WEAPONS_ARR_CNT:NARCOTICS_ARR_CNT | 0.217 | |
| (0.631) | ||
| Constant | 291.547*** | 292.111*** |
| (1.520) | (2.234) | |
| Observations | 9,040 | 9,040 |
| R2 | 0.043 | 0.043 |
| Adjusted R2 | 0.042 | 0.042 |
| Residual Std. Error | 49.688 (df = 9037) | 49.690 (df = 9036) |
| F Statistic | 201.280*** (df = 2; 9037) | 134.213*** (df = 3; 9036) |
| Note: | p<0.1; p<0.05; p<0.01 | |
library(effects)
modEffects <- effect("WEAPONS_ARR_CNT*NARCOTICS_ARR_CNT", intMod)
plot(modEffects)
This is offering the relationship between weapons arrests and score at various levels of narcotics arrests. What we are really looking for are slopes that differ from each other. Knowing that we didn’t find a significant interaction in our model, we probably shouldn’t be surprised that we don’t see any slope difference here.
The five levels that get plotted over is what would lead us towards something like a “floodlight analysis”. Below is more of a “spotlight analysis”:
library(interactions)
interact_plot(intMod, pred = WEAPONS_ARR_CNT, modx = NARCOTICS_ARR_CNT)
crimeScoreGender <-
crimeScore[SEX_CODE_CD != "X", .(SSL_SCORE, SEX_CODE_CD, WEAPONS_ARR_CNT)
][, SEX_CODE_CD := as.factor(SEX_CODE_CD)]
# Below is the dplyr equivalent to what we just used.
# Try it and see how slow it is...smh.
# crimeScoreGender = crimeScore %>%
# filter(SEX_CODE_CD != "X") %>%
# dplyr::select(SSL_SCORE, SEX_CODE_CD, WEAPONS_ARR_CNT) %>%
# mutate(SEX_CODE_CD = as.factor(SEX_CODE_CD))
intMod2 <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT * SEX_CODE_CD, data = crimeScoreGender)
summary(intMod2)
Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT * SEX_CODE_CD, data = crimeScoreGender)
Residuals:
Min 1Q Median 3Q Max
-264.092 -30.093 -0.093 29.908 182.908
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 316.467 7.939 39.862 < 2e-16 ***
WEAPONS_ARR_CNT -2.251 7.202 -0.313 0.75465
SEX_CODE_CDM -16.376 8.005 -2.046 0.04079 *
WEAPONS_ARR_CNT:SEX_CODE_CDM 19.252 7.245 2.657 0.00788 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 53.03 on 19141 degrees of freedom
(379482 observations deleted due to missingness)
Multiple R-squared: 0.02465, Adjusted R-squared: 0.02449
F-statistic: 161.2 on 3 and 19141 DF, p-value: < 2.2e-16Compared to women, men’s score increases by 19.25 on average for each weapons arrest.
Sometimes it helps to see what is going on with a plot:
modEffects <- effect("WEAPONS_ARR_CNT*SEX_CODE_CD", intMod2)
plot(modEffects)
Here is another way (and probably easier) to visualize this effect:
interact_plot(intMod2, pred = WEAPONS_ARR_CNT, modx = SEX_CODE_CD)
The visualization is telling, but the simple slope test is what will tell us which part of the model is significant:
sim_slopes(intMod2, pred = WEAPONS_ARR_CNT, modx = SEX_CODE_CD)
SIMPLE SLOPES ANALYSIS
Slope of WEAPONS_ARR_CNT when SEX_CODE_CD = M:
Est. S.E. t val. p
------- ------ -------- ------
17.00 0.78 21.76 0.00
Slope of WEAPONS_ARR_CNT when SEX_CODE_CD = F:
Est. S.E. t val. p
------- ------ -------- ------
-2.25 7.20 -0.31 0.75We see that the simple slope for F is not significant, whereas the simple slope for M is significant. This offers the “driver” of the relationship.
And another way:
sjPlot::plot_model(intMod2, type = "int",
title = "") +
ggplot2::theme_minimal()
Those plots offered the same information, just slightly different looks – pick whichever one you like most!
For the sake of it, let’s look at one more continuous by continuous interaction:
domNarc <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT * DOMESTIC_ARR_CNT,
data = crimeScore)
summary(domNarc)
Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT * DOMESTIC_ARR_CNT,
data = crimeScore)
Residuals:
Min 1Q Median 3Q Max
-194.551 -30.857 6.657 32.449 224.421
Coefficients:
Estimate Std. Error t value
(Intercept) 276.7127 1.0665 259.462
NARCOTICS_ARR_CNT 2.8284 0.4019 7.038
DOMESTIC_ARR_CNT -1.4681 0.5447 -2.695
NARCOTICS_ARR_CNT:DOMESTIC_ARR_CNT 0.4775 0.2074 2.302
Pr(>|t|)
(Intercept) < 2e-16 ***
NARCOTICS_ARR_CNT 2.02e-12 ***
DOMESTIC_ARR_CNT 0.00704 **
NARCOTICS_ARR_CNT:DOMESTIC_ARR_CNT 0.02132 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 52.34 on 19090 degrees of freedom
(379590 observations deleted due to missingness)
Multiple R-squared: 0.01304, Adjusted R-squared: 0.01288
F-statistic: 84.06 on 3 and 19090 DF, p-value: < 2.2e-16Let’s look at these interactions:
probe_interaction(domNarc, pred = NARCOTICS_ARR_CNT, modx = DOMESTIC_ARR_CNT)
JOHNSON-NEYMAN INTERVAL
When DOMESTIC_ARR_CNT is OUTSIDE the interval [-49.01, -2.40],
the slope of NARCOTICS_ARR_CNT is p < .05.
Note: The range of observed values of DOMESTIC_ARR_CNT is [1.00,
24.00]
SIMPLE SLOPES ANALYSIS
Slope of NARCOTICS_ARR_CNT when DOMESTIC_ARR_CNT = 0.47 (- 1 SD):
Est. S.E. t val. p
------ ------ -------- ------
3.05 0.33 9.33 0.00
Slope of NARCOTICS_ARR_CNT when DOMESTIC_ARR_CNT = 1.60 (Mean):
Est. S.E. t val. p
------ ------ -------- ------
3.59 0.23 15.67 0.00
Slope of NARCOTICS_ARR_CNT when DOMESTIC_ARR_CNT = 2.74 (+ 1 SD):
Est. S.E. t val. p
------ ------ -------- ------
4.14 0.33 12.52 0.00
sjPlot::plot_model(domNarc, type = "int",
title = "") +
theme_minimal()
And also some categorical interactions:
crimeScore2 <- crimeScore[AGE_CURR != "" & SEX_CODE_CD != "X",
AGE_CURR_Rec := .(relevel(as.factor(AGE_CURR), ref = "less than 20"))]
# crimeScore2 = crimeScore %>%
# filter(AGE_CURR != "" & SEX_CODE_CD != "X") %>%
# mutate(AGE_CURR = relevel(as.factor(.$AGE_CURR), ref = "less than 20"))
sexMod <- lm(SSL_SCORE ~ SEX_CODE_CD, data = crimeScore2)
summary(sexMod)
Call:
lm(formula = SSL_SCORE ~ SEX_CODE_CD, data = crimeScore2)
Residuals:
Min 1Q Median 3Q Max
-273.46 -37.69 10.31 42.31 221.31
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 283.4600 0.1868 1517.716 <2e-16 ***
SEX_CODE_CDM -4.7709 0.2145 -22.246 <2e-16 ***
SEX_CODE_CDX -17.2845 7.6793 -2.251 0.0244 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 57.96 on 398681 degrees of freedom
Multiple R-squared: 0.001248, Adjusted R-squared: 0.001243
F-statistic: 249 on 2 and 398681 DF, p-value: < 2.2e-16
ageMod <- lm(SSL_SCORE ~ AGE_CURR_Rec, data = crimeScore2)
summary(ageMod)
Call:
lm(formula = SSL_SCORE ~ AGE_CURR_Rec, data = crimeScore2)
Residuals:
Min 1Q Median 3Q Max
-256.449 -12.741 -1.449 9.567 200.551
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 356.59475 0.09531 3741.5 <2e-16 ***
AGE_CURR_Rec20-30 -35.73399 0.10735 -332.9 <2e-16 ***
AGE_CURR_Rec30-40 -79.14579 0.11226 -705.0 <2e-16 ***
AGE_CURR_Rec40-50 -123.27298 0.12080 -1020.5 <2e-16 ***
AGE_CURR_Rec50-60 -166.16162 0.13211 -1257.7 <2e-16 ***
AGE_CURR_Rec60-70 -207.85339 0.19513 -1065.2 <2e-16 ***
AGE_CURR_Rec70-80 -254.62249 0.47812 -532.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.66 on 398379 degrees of freedom
(298 observations deleted due to missingness)
Multiple R-squared: 0.8958, Adjusted R-squared: 0.8958
F-statistic: 5.708e+05 on 6 and 398379 DF, p-value: < 2.2e-16
sexAge <- lm(SSL_SCORE ~ SEX_CODE_CD * AGE_CURR_Rec, data = crimeScore2)
summary(sexAge)
Call:
lm(formula = SSL_SCORE ~ SEX_CODE_CD * AGE_CURR_Rec, data = crimeScore2)
Residuals:
Min 1Q Median 3Q Max
-256.777 -12.739 -1.384 9.667 200.223
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 349.1487 0.1765 1978.133 < 2e-16
SEX_CODE_CDM 10.4675 0.2093 50.018 < 2e-16
AGE_CURR_Rec20-30 -31.4097 0.2013 -156.055 < 2e-16
AGE_CURR_Rec30-40 -72.7648 0.2144 -339.444 < 2e-16
AGE_CURR_Rec40-50 -116.6520 0.2343 -497.896 < 2e-16
AGE_CURR_Rec50-60 -158.3222 0.2685 -589.698 < 2e-16
AGE_CURR_Rec60-70 -200.0051 0.4588 -435.955 < 2e-16
AGE_CURR_Rec70-80 -246.7251 1.2394 -199.063 < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec20-30 -6.2596 0.2375 -26.356 < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec30-40 -9.0747 0.2513 -36.112 < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec40-50 -9.3964 0.2732 -34.400 < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec50-60 -10.9607 0.3084 -35.540 < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec60-70 -10.9464 0.5072 -21.581 < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec70-80 -10.9950 1.3427 -8.189 2.65e-16
(Intercept) ***
SEX_CODE_CDM ***
AGE_CURR_Rec20-30 ***
AGE_CURR_Rec30-40 ***
AGE_CURR_Rec40-50 ***
AGE_CURR_Rec50-60 ***
AGE_CURR_Rec60-70 ***
AGE_CURR_Rec70-80 ***
SEX_CODE_CDM:AGE_CURR_Rec20-30 ***
SEX_CODE_CDM:AGE_CURR_Rec30-40 ***
SEX_CODE_CDM:AGE_CURR_Rec40-50 ***
SEX_CODE_CDM:AGE_CURR_Rec50-60 ***
SEX_CODE_CDM:AGE_CURR_Rec60-70 ***
SEX_CODE_CDM:AGE_CURR_Rec70-80 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18.56 on 398372 degrees of freedom
(298 observations deleted due to missingness)
Multiple R-squared: 0.8968, Adjusted R-squared: 0.8968
F-statistic: 2.664e+05 on 13 and 398372 DF, p-value: < 2.2e-16
cat_plot(sexAge, pred = SEX_CODE_CD, modx = AGE_CURR_Rec,
geom = "line")
Mediation
As noted earlier, interactions are often referred to as moderators. Just to reiterate, they are testing if the relationship between the X & Y change when introducing a Z. A slightly different model details mediation – the relationship between X and Y exists because M directs the relationship.
We can refer to these paths as follows:
X -> M = a
M -> Y = b
X -> Y = c`
Let’s think how this might look conceptually:
Age -> Salary -> Job Satisfaction
We are saying that the effect of age on job satisfaction goes through salary. In other words, age doesn’t increase job satisfaction – age increases salary, which in turn increases job satisfaction.
To test this mediation model, we need to construct 4 separate regression models:
Model 1: Test path c (X predicts Y)
Model 2: Test path a (X predicts M)
Model 3: Test path b (M predicts Y)
If all of these yield significant results, you can go to:
Model 4: Use X and M as predictor to Y
If the effect of M on Y is significant in Model 4, then there is mediation: full mediation if X is not significant or partial mediation if X is significant. As awesome as a fully-mediated relationship sounds, it is not something that you are going to encounter all that often; instead, partial mediation is where things tend to fall.
Let’s try it out with our crime data:
Just to put some words behind this – we are proposing that narcotics arrests increase weapons arrests, which in turn leads to increased crime scores.
mediationData <- crimeScore[,
.(SSL_SCORE, NARCOTICS_ARR_CNT, WEAPONS_ARR_CNT)]
mediationData <- na.omit(mediationData)
pathC <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT,
data = mediationData)
summary(pathC)
Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT, data = mediationData)
Residuals:
Min 1Q Median 3Q Max
-198.43 -28.12 -2.43 26.57 183.26
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 314.0568 0.8627 364.060 <2e-16 ***
NARCOTICS_ARR_CNT 2.6866 0.2983 9.007 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 50.55 on 9038 degrees of freedom
Multiple R-squared: 0.008896, Adjusted R-squared: 0.008787
F-statistic: 81.13 on 1 and 9038 DF, p-value: < 2.2e-16Our pathC object (Model 1) is significant, so let’s go to the next step:
pathA <- lm(WEAPONS_ARR_CNT ~ NARCOTICS_ARR_CNT,
data = mediationData)
summary(pathA)
Call:
lm(formula = WEAPONS_ARR_CNT ~ NARCOTICS_ARR_CNT, data = mediationData)
Residuals:
Min 1Q Median 3Q Max
-0.2397 -0.2397 -0.2330 -0.2130 4.7670
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.246405 0.008814 141.404 <2e-16 ***
NARCOTICS_ARR_CNT -0.006679 0.003048 -2.191 0.0285 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5165 on 9038 degrees of freedom
Multiple R-squared: 0.000531, Adjusted R-squared: 0.0004204
F-statistic: 4.802 on 1 and 9038 DF, p-value: 0.02846Model 2 – the effect of narcotics on weapons – is also significant. Things are looking promising so far.
pathB <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT,
data = mediationData)
summary(pathB)
Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT, data = mediationData)
Residuals:
Min 1Q Median 3Q Max
-195.052 -28.052 -2.052 26.948 183.948
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 298.215 1.357 219.70 <2e-16 ***
WEAPONS_ARR_CNT 17.837 1.017 17.55 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 49.94 on 9038 degrees of freedom
Multiple R-squared: 0.03294, Adjusted R-squared: 0.03283
F-statistic: 307.8 on 1 and 9038 DF, p-value: < 2.2e-16Model 3 is also looking good! Weapons certainly has an effect on crime score.
pathFull <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT + WEAPONS_ARR_CNT,
data = mediationData)
summary(pathFull)
Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT + WEAPONS_ARR_CNT,
data = mediationData)
Residuals:
Min 1Q Median 3Q Max
-194.221 -28.029 -2.221 26.550 187.586
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 291.5470 1.5197 191.850 <2e-16 ***
NARCOTICS_ARR_CNT 2.8072 0.2933 9.573 <2e-16 ***
WEAPONS_ARR_CNT 18.0598 1.0118 17.849 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 49.69 on 9037 degrees of freedom
Multiple R-squared: 0.04265, Adjusted R-squared: 0.04243
F-statistic: 201.3 on 2 and 9037 DF, p-value: < 2.2e-16Since all predictors continue to be significant, we would have partial mediation. It would be a weak moderation, though, because we did not reduced the magnitude of narcotics when the mediator was included in the model.
We can also test the indirect effect of X through M to Y by subtracting coefficients from models 1 and 4.
# This is Judd and Kenny's method:
pathC$coefficients["NARCOTICS_ARR_CNT"] -
pathFull$coefficients["NARCOTICS_ARR_CNT"]
NARCOTICS_ARR_CNT
-0.1206138 Or by multiplying the full path moderator’s coefficient (model 4) from path a’s coefficient (model 2)
# This is Sobel's method:
pathFull$coefficients["WEAPONS_ARR_CNT"] *
pathA$coefficients["NARCOTICS_ARR_CNT"]
WEAPONS_ARR_CNT
-0.1206138 And those are equal! It doesn’t matter which one you choose, but I find Judd and Kenny to be a little more straightforward.
library(mediation)
medTest <- mediate(pathA, pathFull, sims = 50,
boot = TRUE,
treat = "NARCOTICS_ARR_CNT",
mediator = "WEAPONS_ARR_CNT")
summary(medTest)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.1206 -0.2246 -0.05 <2e-16 ***
ADE 2.8072 2.2102 3.18 <2e-16 ***
Total Effect 2.6866 2.0442 3.09 <2e-16 ***
Prop. Mediated -0.0449 -0.0926 -0.02 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 9040
Simulations: 50 Some of those values should look pretty similar to what we have already done! ACME (Average Causal Mediation Effect) is the indirect effect where we multiplied or substracted our coefficients. ADE (Average Direct Effect) is from our full model. The total effect is what we found from pathC. We are most focused on ACME here, and this significant effect would indicate that we have a significant indirect effect.
We can learn a lot through mediation (it can also be combined with interactions to form mediated moderation or moderated mediation), but it should be used with a great deal of though; it is not just something you start testing models to see how they work. Typically, you need to have a pretty good theoretical reason to dive into mediation.
Effect Sizes
Effect sizes, in conjunction with our p-values, will provide a really good idea about the strength of the difference.
With regard to effect sizes, you will most commonly come across Cohen’s d – it is generally used for t-tests.
Computationally, it is pretty simple:
\[ \frac{\mu_1 - \mu_2}{\sigma}\]
There is also an expanded version:
\[ d = \frac{\mu_1-\mu_2}{\sigma_{pooled}} \]
\[ SD_{pooled} = \sqrt{\frac{\sigma_1^2 + \sigma_2^2}{2}} \]
We are subtracting the mean of one group from another and dividing by the standard deviation.
library(dplyr)
crimeScoreGender[, .(mean = mean(SSL_SCORE),
sd = sd(SSL_SCORE),
n = .N),
by = SEX_CODE_CD]
SEX_CODE_CD mean sd n
1: M 278.6891 59.52397 302320
2: F 283.4600 52.74889 96307
# crimeScoreGender %>%
# group_by(SEX_CODE_CD) %>%
# summarize(mean = mean(SSL_SCORE),
# sd = sd(SSL_SCORE),
# n = n())
sd(crimeScoreGender$SSL_SCORE)
[1] 57.99564We can do it by hand:
(283.46-278.689) / 57.99564
[1] 0.0822648And with the pooled method:
sdPooled <- sqrt((52.74889^2 + 59.52397^2) / 2)
(283.46-278.689) / sdPooled
[1] 0.08483505Or use things already built:
library(compute.es)
tes(t = 23.674, n.1 = 96307, n.2 = 302320)
Mean Differences ES:
d [ 95 %CI] = 0.09 [ 0.08 , 0.09 ]
var(d) = 0
p-value(d) = 0
U3(d) = 53.49 %
CLES(d) = 52.47 %
Cliff's Delta = 0.05
g [ 95 %CI] = 0.09 [ 0.08 , 0.09 ]
var(g) = 0
p-value(g) = 0
U3(g) = 53.49 %
CLES(g) = 52.47 %
Correlation ES:
r [ 95 %CI] = 0.04 [ 0.03 , 0.04 ]
var(r) = 0
p-value(r) = 0
z [ 95 %CI] = 0.04 [ 0.03 , 0.04 ]
var(z) = 0
p-value(z) = 0
Odds Ratio ES:
OR [ 95 %CI] = 1.17 [ 1.16 , 1.19 ]
p-value(OR) = 0
Log OR [ 95 %CI] = 0.16 [ 0.15 , 0.17 ]
var(lOR) = 0
p-value(Log OR) = 0
Other:
NNT = 39.34
Total N = 398627
mes(m.1 = 283.46, m.2 = 278.689,
sd.1 = 52.74889, sd.2 = 59.52397,
n.1 = 96307, n.2 = 302320)
Mean Differences ES:
d [ 95 %CI] = 0.08 [ 0.08 , 0.09 ]
var(d) = 0
p-value(d) = 0
U3(d) = 53.28 %
CLES(d) = 52.32 %
Cliff's Delta = 0.05
g [ 95 %CI] = 0.08 [ 0.08 , 0.09 ]
var(g) = 0
p-value(g) = 0
U3(g) = 53.28 %
CLES(g) = 52.32 %
Correlation ES:
r [ 95 %CI] = 0.04 [ 0.03 , 0.04 ]
var(r) = 0
p-value(r) = 0
z [ 95 %CI] = 0.04 [ 0.03 , 0.04 ]
var(z) = 0
p-value(z) = 0
Odds Ratio ES:
OR [ 95 %CI] = 1.16 [ 1.15 , 1.18 ]
p-value(OR) = 0
Log OR [ 95 %CI] = 0.15 [ 0.14 , 0.16 ]
var(lOR) = 0
p-value(Log OR) = 0
Other:
NNT = 41.96
Total N = 398627Power Analysis
Rules of thumb have been around for a long time and have changed over the years – maybe you learned that you needed 20 rows per predictor, or maybe even 50 rows per predictor. Instead of trusting outdated advice, use actual science to determine how many people you need to find if a difference exists.
We need three of the following parameters:
Effect size
Sample size
Significance level
Power
We should always be doing this a priori.
Power
Power is ability to detect an effect. In NHST words, we are trying to determine if we correctly reject the null hypothesis.
Type I errors: Reject a true \(H_{o}\) (false positive – saying something is there when it is not)
Type II errors: Reject a false \(H_{o}\) (false negative – saying something is not there when it is)
Putting It All Together
Let’s use the pwr package.
library(pwr)
pwr.f2.test(u = 1, v = NULL, f2 = .05, power = .8)
Multiple regression power calculation
u = 1
v = 156.9209
f2 = 0.05
sig.level = 0.05
power = 0.8In the function:
u is the numerator df (k - 1)
v is the denominator df (n - k)
f2 is signficance level
\(\Pi = 1 -\beta\)
\(\beta = Type\,II_{prob}\)
Power is typically set at .8, because it represents a 4 to 1 trade between Type II and Type I errors.
Different Test, Different Power Tests
We just did a test for a linear regression model.
Here is one for a t-test:
tPower <- pwr.t.test(n = NULL, d = 0.1, power = 0.8,
type = "two.sample", alternative = "greater")
plot(tPower)
Let theory be your guide (but be realistic).